// -*- C++ -*-

#ifndef _DG_MATRIX_
#define _DG_MATRIX_

#ifdef DEBUG
#define DG_RANGE_CHECK 1
#endif

#include <assert.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_matrix_double.h>
#include <gsl/gsl_matrix_float.h>
#include <gsl/gsl_vector_float.h>
#include "gsl_linalg_float.h"
#include <omp.h>
 
#include <iostream>
using std::cout;
using std::cerr;
using std::endl;

#include <cmath>
using std::min;

#include <cstring>


// lapack : LU decomposition 

extern "C" void sgesv_(const int&,const int&,float*,const int&,int*,float*,const int&,const int&);
extern "C" void dgesv_(const int&,const int&,double*,const int&,int*,double*,const int&,const int&);

// blas level 1 axpy : y <- a.x + y 
//                     dim        a             x             stride_x   y       stride_y
extern "C" void saxpy_(const int&,const float&, const float*, const int&,float*, const int&);
extern "C" void daxpy_(const int&,const double&,const double*,const int&,double*,const int&);

// blas level 2 gemv : y <- a. op(A) x + b. y

extern "C" void dgemv_( char* tran,const int& m,const int& n,const double& alpha,const double* A,
                        const int& tda,const double* x,const int& incx,const double& b,
                        double* y,const int& incy);

extern "C" void sgemv_( char* tran,const int& m,const int& n,const float& alpha,const float* A,
                        const int& tda,const float* x,const int& incx,const float& b,
                        float* y,const int& incy);

// blas level 3 gemm : C <- c_a op(A).op(B) + c_b C 

extern "C" void sgemm_(char* ta, char* tb,const int& M, const int& N, const int& K, 
                       const float& ALPHA, 
                       const float* A, const int& LDA, 
                       const float* B, const int& LDB,
                       const float& BETA, 
                       float* C, const int& LDC );

extern "C" void dgemm_(char* ta, char* tb,const int& M, const int& N, const int& K, 
                       const double& ALPHA, 
                       const double* A, const int& LDA, 
                       const double* B, const int& LDB,
                       const double& BETA, 
                       double* C, const int& LDC );
  

/*! \brief Flexible data vector format \ingroup dgdatastructure */ 
/*! Features
  <ul>
  <li> proxies
  <li> blas level 1 operations 
  <li> binary operations 
  </ul>
*/

//--------------------------------------------------------------
//VECTOR
//--------------------------------------------------------------
template <
  // types
  class FLOAT,
  class GSL_VECTOR,
  class GSL_VECTOR_VIEW,
  // allocation
  GSL_VECTOR*     (*GSL_VECTOR_CALLOC)       (const size_t),
  void            (*GSL_VECTOR_FREE)         (GSL_VECTOR*),
  // proxies
  FLOAT           (*GSL_VECTOR_GET)          (const GSL_VECTOR*,size_t),
  FLOAT*          (*GSL_VECTOR_PTR)          (GSL_VECTOR*,size_t),
  GSL_VECTOR_VIEW (*GSL_VECTOR_SUBVECTOR)    (GSL_VECTOR*,size_t,size_t),
  // interaction with a scalar
  void            (*GSL_VECTOR_SET)          (GSL_VECTOR*,const size_t,FLOAT),
  void            (*GSL_VECTOR_SET_ALL)      (GSL_VECTOR*,FLOAT),
  int             (*GSL_VECTOR_ADD_CONSTANT) (GSL_VECTOR*,const double),
  int             (*GSL_VECTOR_SCALE)        (GSL_VECTOR*,const double),
  // interaction with another vector
  int             (*GSL_VECTOR_MEMCPY)       (GSL_VECTOR*,const GSL_VECTOR*),
  int             (*GSL_VECTOR_ADD)          (GSL_VECTOR*,const GSL_VECTOR*),
  int             (*GSL_VECTOR_SUB)          (GSL_VECTOR*,const GSL_VECTOR*),
  int             (*GSL_VECTOR_MUL)          (GSL_VECTOR*,const GSL_VECTOR*),
  int             (*GSL_VECTOR_DIV)          (GSL_VECTOR*,const GSL_VECTOR*),
  // blas 1 
  int             (*GSL_BLAS_AXPY)           (FLOAT,const GSL_VECTOR*,GSL_VECTOR*),
  void            (*AXPY)                    (const int&,const FLOAT&, const FLOAT*, const int&,FLOAT*, const int&),
  FLOAT           (*GSL_BLAS_NRM2)           (const GSL_VECTOR*)>
class DGVECTOR
{
  GSL_VECTOR_VIEW view;
public:
  GSL_VECTOR *data;
private:
  size_t r;
public:
  inline size_t size() const {return r;}
  inline size_t size1() const {return r;}
  inline size_t dataTypeSize() const {return sizeof(FLOAT);}
  virtual ~DGVECTOR()  { if(data && data->owner == 1) GSL_VECTOR_FREE (data); }
  DGVECTOR() : data(0) {}
  DGVECTOR(size_t R): r(R) {
    data = GSL_VECTOR_CALLOC (r);
  }
  DGVECTOR (GSL_VECTOR_VIEW _data):view(_data),data(&view.vector),r((_data.vector).size){}
  DGVECTOR(const DGVECTOR&other): r(other.r)
  {
    data = GSL_VECTOR_CALLOC (r);
    GSL_VECTOR_MEMCPY(data, other.data);
  }
  DGVECTOR(size_t n, FLOAT *v) : r(n) {
    data = GSL_VECTOR_CALLOC (r);
    for (unsigned int i=0;i<n;i++)gsl_vector_set(data,i,v[i]);
  }
    
  inline FLOAT operator () (int i) const {
    return GSL_VECTOR_GET (data,i);
  }
  inline FLOAT & operator () (int i) {
    return *GSL_VECTOR_PTR (data,i);
  }
  DGVECTOR & operator = (const DGVECTOR &other)
  {
    if (&other != this)
      {
        if(data)GSL_VECTOR_FREE (data);
        data = GSL_VECTOR_CALLOC (other.data->size);
        GSL_VECTOR_MEMCPY (data, other.data);
        r = other.data->size;
      }
    return *this;
  }

  inline void memcpy(const DGVECTOR& other) 
  {
    GSL_VECTOR_MEMCPY(data,other.data);
  }

  template<class DGVECTOR2>
  inline void memcpy2(const DGVECTOR2& other) 
  {
    for (size_t j=0;j<data->size;j++) data[j*data->stride] = other.data[j*other.data->stride];
  }
  
  inline GSL_VECTOR_VIEW touchSubvector(size_t i,size_t nb)
  {
    return GSL_VECTOR_SUBVECTOR(data,i,nb);
  }
  
  inline const GSL_VECTOR_VIEW seeSubvector(size_t i,size_t nb) const
  {
    return GSL_VECTOR_SUBVECTOR(data,i,nb);
  }
  
  inline void set_all (const double & m ) 
  {
    GSL_VECTOR_SET_ALL (data, m);
  }
  inline void add (const DGVECTOR & y)
  {
    GSL_VECTOR_ADD ( data, y.data );
  }
  inline void sub (const DGVECTOR & y)
  {
    GSL_VECTOR_SUB ( data, y.data );
  }
  inline void add (double d)
  {
    GSL_VECTOR_ADD_CONSTANT (data, d);
  }
  inline void sub (double d)
  {
    GSL_VECTOR_ADD_CONSTANT (data, -d);
  }
  inline void mul (const DGVECTOR & y)
  {
    GSL_VECTOR_MUL ( data, y.data );
  }
  inline void div (const DGVECTOR & y)
  {
    GSL_VECTOR_DIV ( data, y.data );
  }
  inline void scale (const double & y)
  {
    if (y == 0.0)  GSL_VECTOR_SET_ALL( data , 0.);
    else GSL_VECTOR_SCALE ( data, y );
  }
  inline void blas_daxpy (const double & a, const DGVECTOR & y)
  {
    GSL_BLAS_AXPY ( a, y.data, data );
  }
  
  inline void scaleAndAdd(const double& a,const DGVECTOR& y) 
    {
      AXPY(data->size,(FLOAT) a,y.data->data,y.data->stride,data->data,data->stride);
    }
  
  inline double norm(){
    //     return gsl_blas_dnrm2(&view.vector);
    return GSL_BLAS_NRM2(data);
  }
  inline void PRINT (const char *comment = "") const {
    printf ("printing vector %s - size %d\n",comment,(int)size());
    for (size_t i = 0 ; i < size() ; i++){
      printf("%12.5E ",gsl_vector_get (data,i));
      printf("\n");
    }
  }
  
  inline void PRINT_CSTYLE (const char *comment = "") const {
    printf ("%s[%d] = {",comment,(int)size());
    for (size_t i = 0 ; i < size() ; i++){
      if (i!=0) printf(",");
      printf("%f",gsl_vector_get (data,i));
    }
    printf("};\n");
  }
};


typedef DGVECTOR<
  double,gsl_vector,gsl_vector_view,
  gsl_vector_calloc,gsl_vector_free,
  gsl_vector_get,gsl_vector_ptr,gsl_vector_subvector,
  gsl_vector_set,gsl_vector_set_all,gsl_vector_add_constant,gsl_vector_scale,
  gsl_vector_memcpy,gsl_vector_add,gsl_vector_sub,gsl_vector_mul,gsl_vector_div,
  gsl_blas_daxpy,daxpy_,gsl_blas_dnrm2> DGVector;

typedef DGVECTOR<
  float,gsl_vector_float,gsl_vector_float_view,
  gsl_vector_float_calloc,gsl_vector_float_free,
  gsl_vector_float_get,gsl_vector_float_ptr,gsl_vector_float_subvector,
  gsl_vector_float_set,gsl_vector_float_set_all,gsl_vector_float_add_constant,gsl_vector_float_scale,
  gsl_vector_float_memcpy,gsl_vector_float_add,gsl_vector_float_sub,gsl_vector_float_mul,gsl_vector_float_div,
  gsl_blas_saxpy,saxpy_,gsl_blas_snrm2> DGFloatVector;


/*! \brief Flexible matrix structure \ingroup dgdatastructure */ 
/*! Features
  <ul>
  <li> proxies
  <li> blas level 2 and 3 (with DGVector)
  <li> limited lapack support (inversion)
  <li> assignment
  <li> binary operations 
  </ul>
*/

//--------------------------------------------------------------
//MATRIX
//--------------------------------------------------------------
template <
  // types
  class DGVECTOR,
  class FLOAT,
  class GSL_MATRIX,
  class GSL_VECTOR,
  class GSL_MATRIX_VIEW,
  class GSL_VECTOR_VIEW,
  // allocation / deallocation 
  GSL_MATRIX*     (*GSL_MATRIX_CALLOC)      (size_t,size_t),
  void            (*GSL_MATRIX_FREE)        (GSL_MATRIX*),
  // proxies
  FLOAT           (*GSL_MATRIX_GET)         (const GSL_MATRIX*,size_t,size_t),
  FLOAT*          (*GSL_MATRIX_PTR)         (GSL_MATRIX*,size_t,size_t),
  GSL_MATRIX_VIEW (*GSL_MATRIX_SUBMATRIX)   (GSL_MATRIX*,size_t,size_t,size_t,size_t),
  GSL_VECTOR_VIEW (*GSL_MATRIX_COL)         (GSL_MATRIX*,size_t),
  GSL_VECTOR_VIEW (*GSL_MATRIX_ROW)         (GSL_MATRIX*,size_t),
  GSL_VECTOR_VIEW (*GSL_MATRIX_SUBCOL)      (GSL_MATRIX*,size_t,size_t,size_t),
  GSL_VECTOR_VIEW (*GSL_MATRIX_SUBROW)      (GSL_MATRIX*,size_t,size_t,size_t),
  GSL_VECTOR_VIEW (*GSL_MATRIX_DIAGONAL)    (GSL_MATRIX*),
  // operations with scalar 
  void            (*GSL_MATRIX_SET)          (GSL_MATRIX*,size_t,size_t,FLOAT),
  void            (*GSL_MATRIX_SET_ALL)      (GSL_MATRIX*,FLOAT),
  void            (*GSL_MATRIX_SET_IDENTITY) (GSL_MATRIX*),
  int             (*GSL_VECTOR_ADD_CONSTANT) (GSL_VECTOR*,const double),
  int             (*GSL_MATRIX_ADD_CONSTANT) (GSL_MATRIX*,const double),
  int             (*GSL_MATRIX_SCALE)        (GSL_MATRIX*,const double),
  // operations with other matrices
  int             (*GSL_MATRIX_MEMCPY)       (GSL_MATRIX*,const GSL_MATRIX*),
  int             (*GSL_MATRIX_ADD)          (GSL_MATRIX*,const GSL_MATRIX*),
  int             (*GSL_MATRIX_SUB)          (GSL_MATRIX*,const GSL_MATRIX*),
  // blas level 1 
  int             (*GSL_BLAS_AXPY)           (FLOAT,const GSL_VECTOR*,GSL_VECTOR*),
  void            (*AXPY)                    (const int&,const FLOAT&, const FLOAT*, const int&,FLOAT*, const int&),
  // blas level 2 
  int             (*GSL_BLAS_GEMV)           (CBLAS_TRANSPOSE_t,FLOAT,const GSL_MATRIX*,const GSL_VECTOR*,FLOAT,GSL_VECTOR*),
  // blas level 3
  int             (*GSL_BLAS_GEMM)           (CBLAS_TRANSPOSE_t,CBLAS_TRANSPOSE_t,FLOAT,const GSL_MATRIX*,const GSL_MATRIX*,FLOAT,GSL_MATRIX*),
  // linear algebra
  int             (*GSL_LINALG_LU_DECOMP)    (GSL_MATRIX*,gsl_permutation*,int*),
  double          (*GSL_LINALG_LU_DET)       (GSL_MATRIX*,int),
  int             (*GSL_LINALG_LU_INVERT)    (const GSL_MATRIX*,const gsl_permutation*,GSL_MATRIX*),
  // lapack
  void            (*GESV) (const int&,const int&,FLOAT*,const int&,int*,FLOAT*,const int&,const int&),
  void            (*GEMV) (char*,const int&,const int&,const FLOAT&,const FLOAT*,
                           const int&,const FLOAT*,const int&,const FLOAT&,
                           FLOAT*,const int&),
  void            (*GEMM) (char*, char*,const int&, const int&, const int&, 
                           const FLOAT&, 
                           const FLOAT*, const int&, 
                           const FLOAT*, const int&,
                           const FLOAT&, 
                           FLOAT*, const int&)>
            
class DGMATRIX
{
  GSL_MATRIX_VIEW view;
public:
  
  GSL_MATRIX *data;
  size_t* permutation;

  inline size_t size1() const {return data->size1;}
  inline size_t size2() const {return data->size2;}
  inline size_t dataTypeSize() const {return sizeof(FLOAT);}

  virtual ~DGMATRIX() {
    if(data && data->owner == 1)GSL_MATRIX_FREE (data);
    if (permutation) delete [] permutation;
  }

  DGMATRIX(GSL_MATRIX_VIEW _data) : view(_data),data(&view.matrix),permutation(NULL) {}
  
  DGMATRIX(size_t R, size_t C) : permutation(NULL)
  { data = GSL_MATRIX_CALLOC (R, C);  }
  
  DGMATRIX() : data(0),permutation(NULL)
  {}

  DGMATRIX(const DGMATRIX&other):
    data(0),permutation(NULL)
  {
    if(data)GSL_MATRIX_FREE (data);
    data = GSL_MATRIX_CALLOC (other.data->size1, other.data->size2);
    GSL_MATRIX_MEMCPY (data, other.data);
  }

  void memcpy (const DGMATRIX&other)
  {
    GSL_MATRIX_MEMCPY (data, other.data);
  }
  
  template <class DGMATRIX2>
  void copy(const DGMATRIX2 & other)
  {
    for (size_t i=0;i<data->size1;i++) {
      for (size_t j=0;j<data->size2;j++) {
        data->data[i*data->tda+j] = other.data->data[i*other.data->tda+j];
      }
    }
  }
  
  
  DGMATRIX & operator = (const DGMATRIX&other)
  {
    if (&other != this)
      {
        if(data)GSL_MATRIX_FREE (data);
        data = GSL_MATRIX_CALLOC (other.data->size1, other.data->size2);
        GSL_MATRIX_MEMCPY (data, other.data);
      }
    return *this;
  }
  inline FLOAT operator () (int i, int j) const
  {
#ifdef DG_RANGE_CHECK
    return GSL_MATRIX_GET (data,i,j);
#else
    return *((FLOAT *) (data->data + (i * data->tda + j))) ;
#endif
  }
  inline FLOAT & operator () (int i, int j)
  {
#ifdef DG_RANGE_CHECK
    return *GSL_MATRIX_PTR (data,i,j);
#else
    return *((FLOAT *) (data->data + (i * data->tda + j))) ;
#endif
  }
  // data = c_a (x . b) + c_b data
  inline void blas_dgemm(const DGMATRIX & x, const DGMATRIX & b, const FLOAT c_a = 1.0, const FLOAT c_b = 1.0)
  {      
    GSL_BLAS_GEMM (CblasNoTrans,CblasNoTrans, c_a, x.data, b.data, c_b, data);
  }
  
//   openmp parallel version of blas dgemm
//  consider N threads
//   each thread applies ablas_dgemm on a serie of lines
#ifdef MYMY_OPENMP
  inline void blas_dgemm_omp_rows(const DGMATRIX & x, const DGMATRIX & b, const FLOAT c_a = 1.0, const FLOAT c_b = 1.0)
  { 
    int iam = 0, np=1;
#pragma omp parallel private(iam)      
    {
      iam = omp_get_thread_num();
      np  = omp_get_num_threads();
      int istart = size1()*iam/np;
      int iend   = size1()*(iam+1)/np;
      DGMATRIX temp1(touchSubmatrix(istart,iend-istart,0,size2())); 
      DGMATRIX temp2(x.seeSubmatrix(istart,iend-istart,0,x.size2())); 
      temp1.blas_dgemm (temp2,b,c_a,c_b);
    }
  }
#else
  inline void blas_dgemm_omp_rows(const DGMATRIX & x, const DGMATRIX & b, const FLOAT c_a = 1.0, const FLOAT c_b = 1.0)
  {
    GSL_BLAS_GEMM (CblasNoTrans,CblasNoTrans, c_a, x.data, b.data, c_b, data);
  } 
#endif
  
  //   each thread applies ablas_dgemm on a serie of colums
#ifdef MYMY_OPENMP
  inline void blas_dgemm_omp_col(const DGMATRIX & x, const DGMATRIX & b, const FLOAT c_a = 1.0, const FLOAT c_b = 1.0)
  { 
    int iam = 0, np=1;
#pragma omp parallel private(iam)      
    {
      iam = omp_get_thread_num();
      np  = omp_get_num_threads();
      int istart = size2()*iam/np;
      int iend   = size2()*(iam+1)/np;
      DGMATRIX temp1(touchSubmatrix(0,size1(),istart,iend-istart)); 
      DGMATRIX temp2(b.seeSubmatrix(0,b.size1(),istart,iend-istart)); 
      temp1.blas_dgemm (x,temp2,c_a,c_b);
    }
  }
#else
  inline void blas_dgemm_omp_col(const DGMATRIX & x, const DGMATRIX & b, const FLOAT c_a = 1.0, const FLOAT c_b = 1.0)
  {
    GSL_BLAS_GEMM (CblasNoTrans,CblasNoTrans, c_a, x.data, b.data, c_b, data);
  } 
#endif
  
  // mixed data type implementation 
  
//   template <class DGMATRIX2,class DGMATRIX3>
//   inline void gemm(const DGMATRIX2& x,const DGMATRIX3& b,const double c_a = 1.0,const double c_b = 1.) 
//   {
//     size_t tda = data->tda;
//     size_t tdax = x.data->tda;
//     size_t tdab = b.data->tda;
    
//     for (size_t i=0;i<data->size1;i++) {
//       for (size_t j=0;j<data->size2;j++) {
//         FLOAT* dd = data->data + i * tda + j;
//         *dd *= c_b;
//         for (size_t k=0;k<x.data->size2;k++) {
//           *dd += c_a * (x.data->data[i*tdax + k] * b.data->data[k*tdab+j]);
//         }
//       }
//     }
//   }
  
  inline void lu_solve (const DGVECTOR & rhs, DGVECTOR & result)
  {
    int s;
    gsl_permutation * p = gsl_permutation_alloc (size1());
    gsl_linalg_LU_decomp ( data, p, &s);
    gsl_linalg_LU_solve ( data ,  p, rhs.data, result.data ) ;
    gsl_permutation_free (p);
  }

  inline GSL_MATRIX_VIEW column_range (int i, int j)
  {
    return GSL_MATRIX_SUBMATRIX (data,0,data->size1,i,j);
  }

  inline GSL_VECTOR_VIEW touchColumn(int i) 
  {
    return GSL_MATRIX_COL(data,i);
  }

  inline const GSL_VECTOR_VIEW seeColumn(int i) const
  {
    return GSL_MATRIX_COL(data,i);
  }

  inline GSL_VECTOR_VIEW touchRow(int i) 
  {
    return GSL_MATRIX_ROW(data,i);
  }

  inline const GSL_VECTOR_VIEW seeRow(int i) const
  {
    return GSL_MATRIX_ROW(data,i);
  }
  
  inline GSL_VECTOR_VIEW touchSubRow(size_t i,size_t offset,size_t n) 
  {
    return GSL_MATRIX_SUBROW(data,i,offset,n);
  }

  inline const GSL_VECTOR_VIEW seeSubRow(size_t i,size_t offset,size_t n) const
  {
    return GSL_MATRIX_SUBROW(data,i,offset,n);
  }
  
  inline GSL_MATRIX_VIEW touchSubmatrix(int i0,int ni,int j0,int nj) 
  {
    return GSL_MATRIX_SUBMATRIX (data,i0,j0,ni,nj);
  }
  
  inline const GSL_MATRIX_VIEW seeSubmatrix(int i0,int ni,int j0,int nj) const 
  {
    return GSL_MATRIX_SUBMATRIX (data,i0,j0,ni,nj);
  }

  inline void lu_decompose(gsl_permutation*& p,double& det) 
  {
    int s;
    p = gsl_permutation_alloc (size1());
    gsl_linalg_LU_decomp ( data, p, &s);
    det = gsl_linalg_LU_det(data,s);
  }
  
  FLOAT determinant() const 
  {
    DGMATRIX copy = *this;
    FLOAT det;
    copy.invertSecure(det);
    return det;
  }  
  
  inline void invert ()
  {
    int s;
    gsl_permutation * p = gsl_permutation_alloc (size1());
    GSL_LINALG_LU_DECOMP ( data, p, &s);
    GSL_MATRIX *data_inv = GSL_MATRIX_CALLOC (size1(), size2());
    GSL_LINALG_LU_INVERT ( data ,  p, data_inv ) ;
    GSL_MATRIX_MEMCPY (data, data_inv);
    GSL_MATRIX_FREE (data_inv);
    gsl_permutation_free (p);
  }

  inline bool invertSecure (double& det)
  {
    int s;
    gsl_permutation * p = gsl_permutation_alloc (size1());
    GSL_LINALG_LU_DECOMP ( data, p, &s);
    det = GSL_LINALG_LU_DET(data,s);
    GSL_MATRIX *data_inv = GSL_MATRIX_CALLOC (size1(), size2());
    GSL_LINALG_LU_INVERT ( data ,  p, data_inv ) ;
    GSL_MATRIX_MEMCPY (data, data_inv);
    GSL_MATRIX_FREE (data_inv);
    gsl_permutation_free (p);

    return (det != 0.);
  }
  

   bool invert_lapack() {    

     size_t nI = data->size1;
     size_t tda = data->tda;
    
     FLOAT* invA = new FLOAT[nI*nI];
     int * piv = new int[nI];
   
     for (size_t i=0;i<nI*nI;i++) invA[i     ] = 0.;
     for (size_t i=0;i<nI;i++)    invA[i*nI+i] = 1.;
   
     int info;
     
     GESV(nI,nI,data->data,data->tda,piv,invA,nI,info);
    
     if (tda==nI) std::memcpy(data->data,invA,nI*nI*sizeof(FLOAT));
     else for (size_t i=0;i<nI;i++) std::memcpy(data->data+i*tda,invA + i*nI,nI*sizeof(FLOAT));

     delete [] invA;
     delete [] piv;

     return (info == 0) ;
   }

  
  template <class DGVECTOR2>
  inline void mult3(const DGVECTOR2& x, DGVECTOR2& b) const 
  {
    size_t strideB = b.data->stride;
    size_t strideX = x.data->stride;
    size_t tda     = data->tda;
    size_t nI      = b.data->size;
    size_t nJ      = x.data->size;

    // gdot version - should be changed in gaxpy
    // cannot use blas : mixed data types
    
    int iam = 0;
    int np  = 1;
#ifdef MY_OPENMP
#pragma omp parallel private (iam)
    {
      iam = omp_get_thread_num();
      np  = omp_get_num_threads();
      
      size_t iBeg = iam * nI /np;
      size_t iEnd = (iam+1) * nI / np;

#else
      size_t iBeg = 0;
      size_t iEnd = nI;
#endif
      
      for (size_t i=iBeg;i<iEnd;i++) {  
        double *px = x.data->data + i*strideX;
        double *pb = b.data->data + i*strideB;
      
        const FLOAT *pmat = data->data + i*tda;
        for (size_t j=0;j<nJ;j++) *pb -= pmat[j] * *(px++);
      } 
#ifdef MY_OPENMP
    }
#endif
  }  
      
      

  template <class DGVECTOR2>
  inline void mult2(const DGVECTOR2& x, DGVECTOR2& b,const double c_x = 1.,const double c_b=1.0) const 
  {
    size_t strideB = b.data->stride;
    size_t strideX = x.data->stride;
    size_t tda     = data->tda;
    size_t nI      = b.data->size;
    size_t nJ      = x.data->size;

    // gdot version - should be changed in gaxpy
    // cannot use blas : mixed data types
    
    int iam = 0;
    int np  = 1;
#ifdef MY_OPENMP
#pragma omp parallel private (iam)
    {
      iam = omp_get_thread_num();
      np  = omp_get_num_threads();
      
      size_t iBeg = iam * nI /np;
      size_t iEnd = (iam+1) * nI / np;

#else
      size_t iBeg = 0;
      size_t iEnd = nI;
#endif
   
      double *px = x.data->data;
      double *pb = b.data->data;
      
      if (c_b != 1.) for (size_t i=iBeg;i<iEnd;i++) pb[i*strideB] *= c_b;
      
      if (c_x != 0.) {
        if (c_x == 1.) {
            for (size_t i=iBeg;i<iEnd;i++) {
              const FLOAT *pmat = data->data + i*tda;
              for (size_t j=0;j<nJ;j++) pb[i*strideB] += pmat[j] * px[j*strideX];
            }
        }
        else {
          for (size_t i=iBeg;i<iEnd;i++) {
            const FLOAT *pmat = data->data + i*tda;
            for (size_t j=0;j<nJ;j++) pb[i*strideB] += c_x * pmat[j] * px[j*strideX];
          }
        }
      }
#ifdef MY_OPENMP
    }
#endif
  }  
      
  inline void mult (const DGVECTOR & x, DGVECTOR & b, const FLOAT c_x=1.0, const FLOAT c_b=1.0) const
  {// b = c_b * b + c_x * data * x
    GSL_BLAS_GEMV (CblasNoTrans, c_x, data, x.data, c_b, b.data);
  }
  inline void mult (const DGVECTOR& x, DGVECTOR & b, bool transpose,const FLOAT c_x = 1.0,const FLOAT c_b=1.) const 
  {
    CBLAS_TRANSPOSE_t trans = transpose ? CblasTrans:CblasNoTrans;
    GSL_BLAS_GEMV(trans,c_x,data,x.data,c_b,b.data);
  }
  inline void gemv(const DGVECTOR& x, DGVECTOR& b,bool transpose,const FLOAT c_x=1.,const FLOAT c_b=1.) const {
    char trans = transpose ? "n":"t";
    gemv(&trans,data->size1,data->size2,c_x,data->data,data->tda,x->data->data,x->data->stride,c_b,
         b->data->data,b->data->stride);
  }
  inline void add (const DGMATRIX & m ) 
  {
    GSL_MATRIX_ADD (data, m.data);
  }
  inline void add (const double& a ) 
  {
    GSL_MATRIX_ADD_CONSTANT (data, a);
  }
  inline void addDiagonal(double a) 
  {
    DGVECTOR diag(GSL_MATRIX_DIAGONAL(data));
    GSL_VECTOR_ADD_CONSTANT (diag.data,a);
  }
  inline void sub (const DGMATRIX & m ) 
  {
    GSL_MATRIX_SUB (data, m.data);
  }
  inline void scale (const double& m ) 
  {
    if (m == 0.0)  GSL_MATRIX_SET_ALL ( data, 0. );
    else GSL_MATRIX_SCALE (data, m);
  }
  
  inline void scaleAndAdd(FLOAT fact,const DGMATRIX& b) 
  {
    const size_t tda  = data->tda;
    const size_t tda2 = b.data->tda;
    const size_t nI = data->size1;
    const size_t nJ = data->size2;

    if ((tda == nJ) && (tda2 == nJ)) AXPY(nJ*nI,fact,b.data->data,1,data->data,1);
    else {
      for (size_t i=0;i<nI;i++) {
        FLOAT* pa = data->data + i*tda;
        const FLOAT* pb = b.data->data + i*tda2;
        AXPY(nJ,fact,pb,1,pa,1);
      }
    }
  }
  
  template<class DGMatrix2>
  inline void scaleAndAdd_mixed(FLOAT fact,const DGMatrix2& b) {
    
    const size_t tda  = data->tda;
    const size_t tda2 = b.data->tda;
    const size_t nI = data->size1;
    const size_t nJ = data->size2;

    for (size_t i=0;i<nI;i++) {
      for (size_t j=0;j<nJ;j++) {
        data->data[i*tda + j] += fact * b.data->data[i*tda2+j];
      }
    }
  }
    

  
  inline void scaleAndAddTranspose (const double& fact, const DGMATRIX& b) {
#ifdef DG_RANGE_CHECK
    if ((size1() != b.size2()) || (size2() != b.size1())) {
      cerr << "Error in matrix dimensions in function DGMATRIX::scaleAndAddTranspose " << endl;
      exit(1);
    }
#endif
    
    for (unsigned int i=0;i<size1();i++) {
      FLOAT* pa = (data->data)   + i * data->tda;
      FLOAT* pb = (b.data->data) + i;
#ifdef __INTEL_COMPILER
#pragma vector always
#endif
      for (unsigned int j=0;j<size2();j++) {
        *(pa++) += fact * (*(pb));
        pb += b.data->tda;
      }
    }
  }

  inline void scaleAndAddRow(int i,FLOAT fact,const DGVECTOR& vec) 
  {

#ifdef DG_RANGE_CHECK
    if (size2() != vec.size()) {
      cerr << "Error in matrix or vector dimensions in function DGMATRIX::scaleAndAddRow " << endl;
      exit(1);
    }
#endif
    
//     FLOAT *pa = data->data + i * data->tda;
//     FLOAT *pv = vec.data->data;
//     size_t nbCols = size2();
//     size_t stride = vec.data->stride;
//     // for (unsigned int k=0;k<nbCols;k++) pa[k] += pv[k*stride] * fact;
    
    AXPY(data->size2,fact,vec.data->data,vec.data->stride,data->data + i *data->tda,1);
  }

  inline void scaleAndAddCol(int i,FLOAT fact,const DGVECTOR& vec) 
  {
    
#ifdef DG_RANGE_CHECK
    if (size1() != vec.size()) {
      cerr << "Error in matrix or vector dimensions in function DGMATRIX::scaleAndAddCol " << endl;
      exit(1);
    }
#endif
//     FLOAT *pa = data->data + i;
//     FLOAT *pv = vec.data->data;
//     size_t tda = data->tda;
//     size_t stride = vec.data->stride;
    // size_t nbRows = size1();
    // for (unsigned int k=0;k<nbRows;k++) pa[k*tda] += pv[k*stride] * fact;
    AXPY(data->size1,fact,vec.data->data,vec.data->stride,data->data + i,data->tda);
  }

  inline void scaleAndAddColToCol(FLOAT fact,const DGMATRIX& b,int iColB,int iColA) {
    AXPY(data->size1,fact,b.data->data + iColB,b.data->tda,data->data + iColA,data->tda);
  }
  inline void scaleAndAddRowToCol(FLOAT fact,const DGMATRIX& b,int iRowB,int iColA) {
    AXPY(data->size1,fact,b.data->data + iRowB * b.data->tda,1,data->data + iColA,data->tda);
  }
  inline void scaleAndAddColToRow(FLOAT fact,const DGMATRIX& b,int iColB,int iRowA) {
    AXPY(data->size2,fact,b.data->data + iColB,b.data->tda,data->data + iRowA * data->tda,1);
  }
  inline void scaleAndAddRowToRow(FLOAT fact,const DGMATRIX& b,int iRowB,int iRowA) {
    AXPY(data->size2,fact,b.data->data + iRowB * b.data->tda,1,data->data + iRowA * data->tda,1);
  }
  inline void scaleAndAddPartialColToCol(FLOAT fact,int nb,const DGMATRIX& b,int iColB,int iRowB,int iColA,int iRowA) {
    AXPY(nb,fact,b.data->data + iRowB * b.data->tda + iColB,b.data->tda,data->data + iRowA * data->tda + iColA,data->tda);
  }
  inline void scaleAndAddPartialRowToCol(FLOAT fact,int nb,const DGMATRIX& b,int iRowB,int iColB,int iColA,int iRowA) {
    AXPY(nb,fact,b.data->data + iRowB * b.data->tda + iColB,1,data->data + iRowA * data->tda + iColA,data->tda);
  }
  inline void scaleAndAddPartialColToRow(FLOAT fact,int nb,const DGMATRIX& b,int iColB,int iRowB,int iRowA,int iColA) {
    AXPY(nb,fact,b.data->data + iRowB * b.data->tda + iColB,b.data->tda,data->data + iRowA * data->tda + iColA,1);
  }
  inline void scaleAndAddPartialRowToRow(FLOAT fact,int nb,const DGMATRIX& b,int iRowB,int iColB,int iRowA,int iColA) {
    AXPY(nb,fact,b.data->data + iRowB * b.data->tda + iColB,1,data->data + iRowA * data->tda + iColA,1);
  }
  
  inline DGMATRIX& operator+= (const DGMATRIX& b) 
  {
    add(b);
    return *this;
  }
  

  inline DGMATRIX& operator-= (const DGMATRIX& b) {
    sub(b);
    return *this;
  }

  GSL_MATRIX_VIEW submatrix(const int& iRow, const int&nbRow,
                            const int& iCol, const int nbCol)
  {
    return GSL_MATRIX_SUBMATRIX (data, iRow, iCol, nbRow,nbCol);
  }

  const GSL_MATRIX_VIEW seeElementCoefs(const int&,const int&) const;


  inline void set(const double& m,size_t i,size_t j) 
  {
    GSL_MATRIX_SET(data,i,j,m);
  }
  
  inline void set_all (const double& m ) 
  {
    GSL_MATRIX_SET_ALL (data, m);
  }
  inline void set_identity() 
  {
    GSL_MATRIX_SET_IDENTITY(data);
  }

  // change this: expensive !!
  DGMATRIX transpose() const
  {
    DGMATRIX trans(size2(),size1());
    for (unsigned int i=0;i<size1();i++) {
      for (unsigned int j=0;j<size2();j++) {
        trans(j,i) = (*this)(i,j);
      }
    }
    return trans;
  }
  
  DGMATRIX cofactor(int i,int j) const 
  {
    int ni = size1();
    int nj = size2();
    DGMATRIX cof(ni-1,nj-1);
    if (i>0) {
      if (j>0)    DGMATRIX(cof.touchSubmatrix(0,i,0,   j  )).memcpy(DGMATRIX(seeSubmatrix(0,i,0  ,   j)));
      if (j<nj-1) DGMATRIX(cof.touchSubmatrix(0,i,j,nj-j-1)).memcpy(DGMATRIX(seeSubmatrix(0,i,j+1,nj-j-1)));
    }
    
    if (i<ni-1) {  
      if (j<nj-1) DGMATRIX(cof.touchSubmatrix(i,ni-i-1,j,nj-j-1)).memcpy(DGMATRIX(seeSubmatrix(i+1,ni-i-1,j+1,nj-j-1)));
      if (j>0)    DGMATRIX(cof.touchSubmatrix(i,ni-i-1,0,     j)).memcpy(DGMATRIX(seeSubmatrix(i+1,ni-i-1,0  ,   j)));
    } 
    return cof;
  }
  
  void PRINT (const char *comment = "") const
  {
    printf ("printing matrix %s\n",comment);
    for (unsigned int i = 0 ; i < size1() ; i++)
      {
        for (unsigned int j = 0 ; j < size2() ; j++)
          {
            printf("%12.5E ",GSL_MATRIX_GET (data,i,j));
          }
        printf("\n");
      }
    printf("\n");
  }


  inline void PRINT_CSTYLE(const char* comment = "") const {
    printf("%s[%d][%d] = {",comment,(int)size1(),(int)size2());
    for (size_t i=0;i<size1();i++){
      printf("{");
      for (size_t j=0;j<size2();j++) {
        if (j == 0) printf("%1.10f",GSL_MATRIX_GET (data,i,j));
        else        printf(",%1.10f",GSL_MATRIX_GET (data,i,j));
      }
      if (i!=size1()-1) printf("},\n");
      else              printf("}\n");
    }
    printf("};\n");
  }
};

typedef DGMATRIX<DGVector,double,gsl_matrix,gsl_vector,gsl_matrix_view,gsl_vector_view,
                 gsl_matrix_calloc,gsl_matrix_free,
                 gsl_matrix_get,gsl_matrix_ptr,
                 gsl_matrix_submatrix,
                 gsl_matrix_column,gsl_matrix_row,gsl_matrix_subcolumn,gsl_matrix_subrow,
                 gsl_matrix_diagonal,
                 gsl_matrix_set,gsl_matrix_set_all,gsl_matrix_set_identity,
                 gsl_vector_add_constant,gsl_matrix_add_constant,gsl_matrix_scale,
                 gsl_matrix_memcpy,gsl_matrix_add,gsl_matrix_sub,
                 gsl_blas_daxpy,daxpy_,gsl_blas_dgemv,gsl_blas_dgemm,
                 gsl_linalg_LU_decomp,gsl_linalg_LU_det,gsl_linalg_LU_invert,dgesv_,dgemv_,dgemm_> DGMatrix;

typedef DGMATRIX<DGFloatVector,float,gsl_matrix_float,gsl_vector_float,
                 gsl_matrix_float_view,gsl_vector_float_view,
                 gsl_matrix_float_calloc,gsl_matrix_float_free,
                 gsl_matrix_float_get,gsl_matrix_float_ptr,
                 gsl_matrix_float_submatrix,
                 gsl_matrix_float_column,gsl_matrix_float_row,gsl_matrix_float_subcolumn,gsl_matrix_float_subrow,
                 gsl_matrix_float_diagonal,
                 gsl_matrix_float_set,gsl_matrix_float_set_all,gsl_matrix_float_set_identity,
                 gsl_vector_float_add_constant,gsl_matrix_float_add_constant,gsl_matrix_float_scale,
                 gsl_matrix_float_memcpy,gsl_matrix_float_add,gsl_matrix_float_sub,
                 gsl_blas_saxpy,saxpy_,gsl_blas_sgemv,gsl_blas_sgemm,
                 gsl_linalg_LU_decomp_float,gsl_linalg_LU_det_float,gsl_linalg_LU_invert_float,sgesv_,sgemv_,sgemm_> DGFloatMatrix;


/*! \brief Outer product of two DGVectors resulting in a DGMatrix \ingroup dgdatastructure */ 
// should be rather slow, maybe we can access the elements directly 

template <class DGMATRIX,class DGVECTOR>
inline void outer_prod ( const DGVECTOR & v1  , const DGVECTOR & v2, DGMATRIX & m )
{
  for (size_t i=0;i<v1.size();i++) {  
    for (size_t j=0;j<v2.size();j++) {
      m.set(v1(i)*v2(j),i,j);
    }
  }
  
  //   for (unsigned int i=0;i<v1.size();i++)
  //     for (unsigned int j=0;j<v2.size();j++)
  //       gsl_matrix_set (m.data,i,j,v1(i)*v2(j));
}

class DG3x3Matrix : public DGMatrix
{
public:
  DG3x3Matrix () : DGMatrix(3,3) { }  
};

#endif
